Resonant forcing of nonlinear systems of differential equations 
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We study resonances of nonlinear systems of differential equations, including but not limited to 
the equations of motion of a particle moving in a potential. We use the calculus of variations to 
determine the minimal additive forcing function that induces a desired terminal response, such as 
an energy in the case of a physical system. We include the additional constraint that only select 
degrees of freedom be forced, corresponding to a very general class of problems in which not all of 
the degrees of freedom in an experimental system are accessible to forcing. We find that certain 
Lagrange multipliers take on a fundamental physical role as the effective forcing experienced by 
the degrees of freedom which are not forced directly. Furthermore, we find that the product of the 
displacement of nearby trajectories and the effective total forcing function is a conserved quantity. 
We demonstrate the efficacy of this methodology with several examples. 

PACS numbers: 05.45.Xt, 05.45.-a 
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Resonance in nonlinear systems is an impor- 
tant topic that has been explored in depth [1, 
[H, 0|. Resonance in a linear system is de- 
fined to be a maximum response amplitude when 
driven by a signal at a specific frequency. In 
this case something about the forcing function, 
namely, the forcing frequency, mirrors something 
about the system, its natural frequency. Previ- 
ous work on driven damped nonlinear oscillators 
demonstrated that such a system will achieve a 
maximum amplitude when the forcing dynamics 
matches the time-reversed dynamics of the same 
system without forcing [4]. Again there is a rela- 
tionship between the natural dynamics of the sys- 
tem and the dynamics of the drive. In this paper 
we derive the most efficient resonant forcing func- 
tion possible for a very general class of systems, 
namely, systems which can be described by cou- 
pled first-order differential equations, including 
systems which exhibit chaos [5]. In the method- 
ology we present, there is no restriction on the 
degrees of freedom which may be forced, so it is 
possible to compute the resonant forcing of, say, 
one of two coupled oscillators. In such a system, 
only one of the four degrees of freedom would 
be forced; this was not possible previously. We 
show that optimal forcing functions may be used 
for system parameter identification via resonance 
spectroscopy. Furthermore, conservation laws in 
closed systems usually correspond to a fundamen- 
tal symmetry. In this paper we show that an 
open dissipative system subject to optimal reso- 
nant forcing has a special conserved quantity and 
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a corresponding symmetry. This conserved quan- 
tity is the dot product of the separation of nearby 
trajectories and the effective forcing experience 
by all degrees of freedom. 



I. INTRODUCTION 

There has been extensive work on sinusoidally driven 
nonlinear oscillators in the contexts of synchroniza- 
tion stochastic resonance @, H[ and nonlinear re- 
sponse phenomena |9|, |10| . Resonance phenomena of non- 
linear systems due to aperiodic and chaotic forcing func- 
tions [3, [ll| has been less studied, but results from work 
in this area indicate that generally a nonlinear oscilla- 
tor will have a greater response when driven with the 
correct aperiodic signal rather than a sinusoidal one. A 
related topic is system identification via resonance curves 
of nonlinear systems [l2| and periodically driven chaotic 
systems [La ]. Plapp and Hiibler [l4j and others [l[ have 
used the calculus of variations to show that a special 
class of aperiodic driving forces can achieve a large en- 
ergy transfer to a nonlinear oscillator. Such nonsinu- 
soidal resonant forcing functions yield a high signal-to- 
noise ratio which can be used for high-resolution system 
identification [l^|. In a recent paper, Gintautas, Foster, 
and Hiibler f|J explored resonant forcing of time-discrete 
chaotic dynamics. In this work, we extend their method 
to time- continuous systems of ordinary differential equa- 
tions, including but not limited to, the equations of mo- 
tion of a particle in a potential, and show that the op- 
timal forcing function induces a desired response more 
efficiently than a sinusoidal forcing function. 

Systems of first order differential equations are ubiq- 
uitous in modern science and engineering. Furthermore, 
any higher order differential equation, such as an equa- 
tion of motion for a Hamiltonian system, or system of 
equations may be cast as a set of first order equations. 
Systems of differential equations have been used to model 
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a rich variety of systems, ranging from complex net- 
works [l6| to jet flow |17| , to give two very recent ex- 
amples. In these cases the correct model accurately re- 
produces the natural dynamics of the system. In other 
cases, it is important not only to correctly model the un- 
perturbed dynamics of the system but also to be able to 
control or influence these dynamics. 

In this paper, we present a methodology for determin- 
ing the resonant forcing of a system of first order differ- 
ential equations in which only select degrees of freedom 
are forced. This is motivated by the difficulty or impos- 
sibility of forcing all of the degrees of freedom in certain 
experiments. For example, consider a physical oscillator 
in which it is possible to directly force the position x but 
not the velocity x. Therefore, the method we present 
may be applied to a very general class of problems. We 
show analytically that the resonant forcing functions are 
closely related to the unperturbed dynamics of the sys- 
tem in that the product of the displacement of nearby 
trajectories and the effective total forcing function is a 
conserved quantity. We also show that the optimal forc- 
ing for a damped oscillator moving in a potential is pro- 
portional to the time reflected dynamics of the corre- 
sponding unperturbed system; this is the "principle of 
the dynamical key" explored by Hiibler et al. [2, [4]. Fur- 
thermore, we find that certain Lagrange multipliers take 
on a fundamental physical role as the efficiency of the 
forcing function and the effective forcing experienced by 
the degrees of freedom which are not forced directly. We 
demonstrate the efficacy of the methodology with sev- 
eral examples. Since the method we present is general 
and requires only access to one degree of freedom, nearly 
any system that is accurately modeled using a system of 
first order equations can also in principle be controlled 
efficiently, including systems which exhibit chaos [B[ . 



II. GENERAL FORMULATION 

We begin with a multidimensional first order system 
with forcing: 



x = /(a 



F, 



(1) 



where 



x{t) e 



denotes the state of the d- 



dimensional system at time t, and F = F(t) G ~R d de- 
notes the forcing function at time t. This system has d 
degrees of freedom. We seek to minimize the total forcing 
effort, that is, the integral of the magnitude of F from 
t = to t = t, which we define to be the constant F 2 : 



[F(t) ■ F(t)]dt. 



(2) 



Here the terminal time r is a free parameter. We require 
that < d u < d degrees of freedom be unforced. With- 
out loss of generality, we choose to order the variables 
so that xi, ■ ■ ■ , Xd u are unforced and Xd u +i, ■ ■ ■ , Xd are 



forced. Thus we will require that 

Fi(t) = 0, for i = 1, . . . , du and < t < t, (3) 

where Fi(t) is the ith component of F(t). This prob- 
lem can be solved by a variation of the functional S = 
J T L g dt. The La grange function L g is given by 

L g = L{x,S,F,t) + \K(x,x',t)5 D (t ~ r), (4) 

where <$o(i — r) is the Dirac delta function and A is a 
constant Lagrange multiplier. The function K is a gen- 
eralized boundary condition for t = r and represents a 
constraint at the terminal time. We will require that K 
be in the form 

K[x(t),x*(t),t] = att = r. (5) 

In light of the above constraints, L is given by 

L=±F-F + F-f + jl(t) ■ [Z-f(x)-F]. (6) 

Because the equation of motion in Eq. ([T]) is a nonintegral 
constraint, fl(t) is a time dependent Lagrange multiplier. 
We have defined the vector T(t) = 7i(*)%) where 

7i(i), . . . , 7d„ (t) are time dependent Lagrange multipliers 
and ej is the unit basis vector in the direction of Xj. 
This term represents the constraint that certain degrees 
of freedom not be forced. Thus the Lagrange problem is 



SS = 6 [ L + XK6 D (t-T)dt, 
Jo 



(7) 



where SS is the variation of S. Following Wargitsch and 
Hiibler p}, we derive the Euler-Lagrange equations for 
this problem (in which the terminal time is a free param- 
eter) in the Appendix. The equations of motion are 



dL d f dL\ 
dxi dt V dii ) 
dL 



dF 



0. 



for i = 1, . . . , d. At the upper boundary, for t = r, 

dx; 



BK dL n 

OXi OXi 

^dK T A/. dL\ n 

/=i 



(8) 
(9) 

(10) 
(11) 

(12) 



for i = 1, . . . , d. At the lower boundary, we have the 
initial condition x(0). The equations of motion yield: 



J T /i + /j= 0, 
F + f - p = 0, 



(13) 
(14) 
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where L. 



(dfi/d 



ated at x(t). The superscript T indicates the transpose 
operator. We now define the quantity 



is the Jacobi matrix evalu- 



G = F 



(15) 



noting that the components of G in the direction of un- 
forced degrees of freedom are the corresponding compo- 
nents of r and the components of G in the direction of 
forced degrees of freedom are equal to the correspond- 
ing components of F. When we solve for the Lagrange 
multiplier ft using Eq. (fl4|) , 



then Eqs. © and 



P(t) - G(t), 
reduce to simply 



G 



-3 T G. 



(16) 



(17) 



From Eq. (fTS"]) and (jTTJ) we identify G as the effective 
total forcing function; it reduces to the optimal forcing 
F when we remove the constraint in Eq. ([3]). We iden- 
tify the Lagrange multipliers 71 , . . . 7d„ to be the effective 
forcing experienced by the degrees of freedom j for which 
Fj = 0; this changes the trajectories of these degrees of 
freedom via the coupling in f(x) rather than direct ad- 
ditive forcing via F . We now further simplify the upper 
boundary conditions in Eqs. (jTTJ) and (|12p . which become 



dK 



(18) 



A if = ¥ {t) • P(t) + 6{t) • ^ (r) ] ' (i9) 

for i = 1, . . . , d. These boundary conditions along with 
the initial condition x(0) and Eqs . fTI and [T7l form a com- 
plete boundary value problem which, in principle, may 
be solved analytically or numerically to determine T(t) 
and F(t) for < t < r. 

The control is stable if, on average, the displacement 
of nearby trajectories decreases. Consider a trajectory 
given by Eq. ([T]), and a nearby trajectory given by x! = 

f(x') + F, where x and x 1 are related by e = x — x! . If 
we Taylor expand f(x) for small e, we obtain 



e = Je. 



(20) 



Multiplying both sides of the transpose of Eq. (JTTJ) by e, 
we have G T e — — G T Je. Using Eq. (120]) . this becomes 



G-e 



-G ■ e, or 



dt 



G) =0, 



(21) 



a quantity that is invariant for all t. We define this to be 
the conserved quantity P: 



P = e-G, 



(22) 



and note that P depends on the observables x and F as 
well as the Lagrange multipliers in F, which we have iden- 
tified as the effective indirect forcing of certain degrees 
of freedom. This further reinforces the idea that G rep- 
resents the effective forcing experienced by the system, 
taking into account the coupling via f(x). Note that P 
is conserved even if the unperturbed dynamics is chaotic 
or periodic. We can independently show that P is a con- 
served quantity using the invariance of the Lagrangian. 
Consider the transformation, 



x — ► x + e, 
x — > x + e. 



(23) 
(24) 



Under this transformation, the variation of the La- 
grangian is in the form L — > L + 8L, with 

SL = ft- (t- Je). (25) 

Using Eq. (|20|) . SL = and we immediately see that 
the Lagrangian is invariant under this transformation. 
Noether's theorem [18| states that if 8L for a given trans- 
formation can be written as a total derivative of a func- 
tion U, that is, SL = dll/dt, then there is a corresponding 
conserved quantity j (called Noether's current): 



J 



d 4-u. 

dx 



(26) 



In this case, SL = dU/dt = 0, so U is some constant c, 
and it follows that 



dt 3 dt 



dL 

— - c ) = 0. 

dx 



(27) 



From Eqs. ([6]) and (|T4]) . we have dL/dx = ft = G and 
recover dP/dt = 0. Conservation laws in closed systems 
usually correspond to a fundamental symmetry. Here we 
have shown that an open dissipative system subject to 
optimal resonant forcing has a special conserved quantity 
and a corresponding symmetry. 



III. EXAMPLES 
A. One-dimensional damped oscillator 

We now illustrate the methodology with several ex- 
amples. Consider first a one-dimensional damped driven 
oscillator, 



dV 

x + r ] x + — = F(t), 
ox 



(28) 



where r/ is the coefficient of linear damping and V(x) is a 
time-independent potential. The energy of the oscillator, 
E(t) = x(t) 2 /2 + V[x(t)], will provide a constraint at 
t = t. We first consider a system of two coupled first- 
order equations: 



X-2 



-r\x 2 



av 
dx i 



(29) 
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with initial conditions 



xi(0) 
x 2 (0) 



xo, 



(30) 
(31) 



and require that only F 2 be forced, that is, = 

for all t. We will first solve for the forcing function in 
terms of x(t), then explore the meaning of the conserved 
quantity for this type of system. 



1. Equations of motion and general solution 

Accordingly, T(t) — j(t)ei, where we have defined 
j(t) = 71 (i). This system is equivalent to Eq. (|2"8"j) if 
we identify F 2 (t) = F(t). The Jacobi matrix for this 
system is 




(32) 



For the upper boundary condition, we will use 



K[x(t),S(t),r] = ±a$(t) + V[x 1 (t)]-E = 0, (33) 

where E is a constant energy value we wish the oscillator 
to attain at t = r. Eq. Q17p gives equations of motion for 
F 2 (t) and 7 (t): 



d 2 V 

F 2 (t)= V F 2 (t)-~/(t). 



At the upper boundary, evaluating Eqs. (ITS]) and l[T9l 

gives 



(34) 
(35) 



dV 



F 2 (r) 



t = T 

Xx 2 (t), 



-F 2 (t) +1 x 2 (t)=F 2 (t) 



dV 



dxi 



rjx 2 (r) 



(36) 
(37) 
(38) 



Eqs. (f36jl . (|37| . and ([38| can be solved for A, 7(r), and 
F 2 (t) in terms of xi(r) and £2(1") to provide explicit 
upper boundary conditions: 



A = -277, 
F 2 (r) = 277x2 (r), 

ay 

7(T) = 



(39) 
(40) 

(41) 



Eqs. (E3> (EH), I0U, 601), and gT]) form a well- 

posed boundary value problem. Eliminating 7 from 
Eqs. (|3"4"| and l|33|) gives an equation of motion for F 2 : 



F 2 {t) - V F 2 (t) + F 2 (t) 



d 2 V 
dx\ 



0. 



(42) 



A trial solution for Eqs. ([28)) and (|34|) is given in the form 

F 2 {t)=ar)X 2 (t). (43) 

Using Eq. (|4T)|) , we find that this is a valid solution only 
for a — 2. This is the same result as calculated by War- 
gitsch and Hubler [lj using a different formulation. 



2. Conserved quantity for one- dimensional damped 
oscillator 



Using this solution for F(t), we consider the conserved 
quantity P. Using Eq. (|20p. we obtain the following equa- 
tion of motion for e: 



d 2 V 



(44) 



where we have defined e = ei and eliminated e 2 = t\. 
This equation of motion is valid on the domain < t < 
t. We substitute Eq. (|4*3"|) into Eq. and operate on 
the resulting equation with an additional time derivative 
(henceforth for this example we will use x\ — > x and 
x 2 — > x): 



x(t) - f]x(t) + x(t) 



d 2 V 



dx 2 



= 0. 



(45) 



Under the transformation t —> r—t (from which it follows 
that d/dt — > —d/dt), this equation becomes 



'x(t — t) + rjx(r — t) + x(t — t) 



d 2 v 



dx 2 



= 0, (46) 



t = T-t 

where the second partial derivative d 2 V/dx 2 is evaluated 
at t = T—t. This equation is precisely in the same form as 
Eq. the equation of motion for e. Furthermore, it is 
valid on the same domain, namely, < t < r. Therefore 
we identify 



e(i) = Ax(t - t) oc F 2 (t -t), 



(47) 



with A an arbitrary constant. We use Eqs. (| 1 7[) . (I22p . 
I, and (|4^)) to write P in terms of x: 



P = r)x(t)x(T-t)+x(t) 



dV 

dx 



-x(r-t) 



8V_ 

dx 



, (48) 



which is completely symmetric under the transformation 
t — > r — t. We have absorbed any multiplicative constants 
into P. Since P does not change in time, Eq. (|48|) must 
hold for i = r/2: 



P = r#(|) i 



d 
df 



1 



dt 



(49) 



Here dE(t)/dt is the instantaneous rate of energy energy 
change of the oscillator and is equal to the instantaneous 
rate of energy transfer of the force F 2 (i) . When evaluated 
at t = r/2, it is equal to the conserved quantity P. 



5 



B. Explicit example: isotonic harmonic oscillator 

We now illustrate the above formulation with sev- 
eral explicit examples. Consider a forced isotonic har- 
monic oscillator of the form of Eq. (|'28p. with a potential 
V = x 2 uj 2 /2 + k/x 2 . Thus the equation of motion of the 
oscillator is 



0.40 



x + rji + lo 2 x t = F(t) 



(50) 



We will use generalized initial conditions, with x(0) = .to 
and x(0) = Vq. This potential represents a harmonic 
oscillator with a centripetal barrier [lj| [2(| and the cor- 
responding equation of motion for r\ — and no forcing 
is a particular case of the Pinney-Ermakov equation [2l[ . 
The unforced rj = case is is an example of a nonlinear 
isochronous system, that is, the amplitude of the oscil- 
lations of the solution are independent of the frequency. 
For these examples, however, we will consider the system 
with damping so that rj is left as a free parameter. First, 
we examine the case where k — 0, corresponding to a 
simple damped harmonic oscillator. For optimal forcing, 
F(t) = 2r]x(t), and the solution for x(t) is 



x(t) 



oVt/2 



X() 



cosh | \/ r\ 2 — Au) 2 



(2v - rjxp) 
\Jr\ 2 - 4w 2 

It follows that 

Fit) = 2r)i(t) = 2r 1 e qt/2 



sinh | \Jr\ 2 — 4w 2 



(51) 



cosh I \Jif - Auj 2 



[r)v - 2lo 2 xq) 
yfrj 2 - Auj 2 

The solution for e(t) is 



sinh | \/ r/ 2 — Auj 2 



(52) 



=(*) = e(0)e-" t / 2 cosh § y/ V 2 - 4w 2 



+ [2e(0) + rye(0)] sinh \ ^rj 2 - 4w 2 , (53) 

where the initial conditions for e at t = are e(0) and 
e(0). The conserved quantity can be computed exactly 
to be 




1.99 2.00 2.01 2.02 2.03 2.04 



FIG. 1: Ratio of total forcing effort for optimal forcing to 
that of sinusoidal forcing. Since the optimal forcing is more 
efficient, this ratio is always less than 1 (dashed line). Here, 
z(0) = 1.0, x(0) = 0.01, r) = 0.3, and E = 2.0. A different 
value of r was calculate for each value of ui plotted. 



t = t using less overall effort. When we plot the ratio in 
Fig. Q] for a range of w, we see that it is always less than 
1, as expected. 

At this point we remove the restriction that k = and 
consider the system with a nonperturbative nonlinear- 
ity (k oc w). Since optimal forcing functions for maps 
have been demonstrated to be useful for resonance spec- 
troscopy [2| , we will show this to be the case for the forc- 
ing functions we have calculated above. For a harmonic 
oscillator (with k = 0) , the optimal forcing was calculated 
as a function of the natural frequency u>. Then a test sys- 
tem harmonic oscillator with lo — u>o was forced with this 
function for different values of u> until the energy of the 
oscillator reached the desired value E. The forcing effort 
Fit) ■ Fit) was then integrated from t = to t = r to 
obtain the total forcing effort F 2 . As expected, the ratio 
of response to total forcing effort is maximal when the 
natural frequency of the forcing function matches that of 
the test system [see Fig. |2Ja)]. We repeat this analysis 
with an isotonic harmonic oscillator with a nonlinearity, 
that is, k 7^ 0, and observe similar results [see Fig. Hfb)]. 
Thereby the optimal forcing in terms of an unknown and 
variable parameter may be used for system identification. 



P = uj 2 x(0)e(0)+x(0)i{0). (54) 

We can compare the effectiveness of this forcing function 
to that of sinusoidal driving. For a given ui, we compute 
the time r such that the energy of the oscillator under 
optimal forcing [see Eq. ([52]) ] is equal to the desired value 
given in Eq. (|33|) . Then we force the same oscillator using 
instead F S inusoidai(i) = Asinwt, and choose A such that 
the energy of the system reaches the same value at t = r. 
Then using Eq. ([2|) we compare F 2 for both optimal and 
sinusoidal forcing. We expect that the optimal forcing is 
able to cause the system to reach the desired energy at 



C. Explicit example: linear ODE system 

The method presented may also be applied to non- 
Hamiltonian systems of ordinary differential equations 
(ODEs). Since any ODE system can be written as an 
equivalent hrst-order system, the method is very general. 
We illustrate this with a simple example. We consider the 
following system: 

f ±i\ (axi + kx 2 \ . {FA f ..s 
{x 2 ) = {kx 1+ ax 2 ) + {F 2 )- ^ 
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in terms of the initial conditions, 



1.8735 
1 



2.02 



FIG. 2: Resonance curve for isotonic harmonic oscillator. 
As expected, the ratio of response to total forcing effort is 
maximal when the natural frequency of the forcing function 
matches that of the test system (uj = luq). Here, x(0) = 1.0, 
x(0) = 0.01, i) = 0.3, loo = 2.00, and E - 3.6. (a) Simple har- 
monic oscillator with k — 0. (b) Isotonic harmonic oscillator 
with k — 1. 



and require that only x 2 be forced, that is, d u = 1. The 
Jacobi matrix of this system is constant and symmetric: 



a k 
k a 



(56) 



The method gives rise to the following equations of mo- 
tion: 



7 = —aj - kF 2 , 
F 2 = -k-f - aF 2 . 



(57) 
(58) 



Consider a system in which one degree of freedom, x 2 , 
is accessible to forcing but is coupled to another degree 
of freedom, x 2 , over which we have no control. x\ may 
represent, say, the contact age degree of freedom in a 
sliding friction model [22| . Suppose we want x 2 to reach 
some desired value at t = r but have no control over X\. 
In such a case, at t = r we have the simple boundary 
condition, 



K[x(t),x(t),t] = x 2 (t) 



C p = 0. 



(59) 



We will use this boundary condition for this simple exam- 
ple of two coupled linear ordinary differential equations. 
Then Eqs. (fT8f and (jT9|) give rise to the following explicit 
boundary conditions on 7 and F 2 , as well as the explicit 
value of A: 

= 0, (60) 
F 2 (t) = -2[ax 2 (r) + feci(r)] (61) 
A = -F 2 (t) = 2 [aC p + kxi (t)] (62) 

It is possible to solve the corresponding boundary value 
problem analytically. We write the explicit form of F 2 (t) 



F 2 (t) = -2e~ at (sechfci) coshfc(£ - r) 
x I \kxi(Q) + ax 2 (0)] coshfci 

+ [axi(0) + kx 2 (0)] sinhfct j, 

and we also find an explicit expression for F 2 , 



(63) 



F 



{[feEi(0) + ax 2 (0)] coshfci 



(sech fct) 4 
2(a 3 - afc 2 ) 

[oaii(O) + kx 2 (0)] smhktl 



-k 2 - 



'(2a 2 



k 2 ) 

a 2 cosh 2kt — ak sinh 2kt | . 



(64) 



As with the isotonic harmonic oscillator, we may 
use the calculated forcing function for resonance spec- 
troscopy. For this system the optimal forcing was calcu- 
lated as a function of the coupling parameter k. Then 
a test system with k = ko was forced with this function 
for different values of k until x 2 reached the value C p . 
The forcing effort Fit) ■ F(t) was then integrated from 
t = to t = t to obtain the total forcing effort F 2 . 
As expected, the ratio of response to total forcing effort 
is maximal when the coupling parameter of the forcing 
function matches that of the test system (see Fig. [3]). 
Just as before, the optimal forcing in terms of an un- 
known and variable system parameter can be used for 
system identification. 



0.080 



0.075 
^ 0.070 



0.065 



0.060 




FIG. 3: Resonance curve for linear system of first order differ- 
ential equations. As expected, the ratio of response to total 
forcing effort is maximal when the coupling parameter of the 
forcing function matches that of the test system (k = fco)- 
Here, xi(0) = 1.0, x 2 (0) = 4.0, a = 1.3, k = 0.2, and 
C p = 1.5. 



IV. CONCLUSIONS 

We study resonances of forced systems of ordinary 
differential equations. We use a constraint at termi- 
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nal time [Eq. ([5])] and seek the forcing function which 
minimizes the total effort [Eq. (J2J)], subject to the ad- 
ditional constraint that certain degrees of freedom are 
not directly forced [Eq. ©]. To determine this forcing 
function, we seek the stationary points of the Lagrange 
function [Eq. (|6|)] and thereby obtain equations which de- 
termine the dynamics of the forcing function [Eqs. I|17p- 
(|19p]. From these equations we identify the effective total 
forcing to be a vector comprising the direct forcing and 
the Lagrange multipliers that represent the effective in- 
direct forcing of certain degrees of freedom [Eq. (fT5|) ] . 
We demonstrate that the product of the effective forcing 
and the displacement of nearby trajectories is a conserved 
quantity [Eq. ([22])] . The methodology presented can be 
applied to a very general class of problems in which not 
all of the degrees of freedom in an experimental system 
are accessible to forcing. Furthermore, the methodology 
is not restricted to Hamiltonian systems or systems with 
small forcing but can applied to any system of ordinary 
differential equations. 



APPENDIX: VARIATIONAL PRINCIPLE WITH 
FREE TERMINAL TIME 



Here we derive the equations of motion for a variational 
problem of the form given in Eq. ([7]): 



SS = S L(x\x\t) + \K(x\x\t)6 D (t-T)dt, (A.l) 
Jo 



where x l represent generalized coordinates, with i — 
1, . . . , N. Since the terminal time r is not fixed but is 
a free parameter, we use a parametric representation of 
the problem. Accordingly we replace i, x 1 , and x 1 with 
the following substitution rules: 



We demonstrate the effectiveness of the methodology 
with several examples. We compare forcing calculated 
using a variational principle to sinusoidal forcing for a 
damped harmonic oscillator and find that the sinusoidal 
forcing is less efficient (see Fig. [T]) . We present a reso- 
nance curve for a damped harmonic oscillator as well as a 
nonlinear isotonic harmonic oscillator in Fig. [5] and verify 
explicitly that the optimal effective forcing complements 
the separation of nearby trajectories [Eq. (|54|) ]. We also 
apply this method to a forced linear system of first order 
differential equations [Eq. ([55)) ]. We solve for the exact 
optimal forcing as a function of the terminal time ana- 
lytically demonstrate that the solution gives the correct 
peak in the resonance curve (see Fig. [3]) . Thus we show 
that the method may be used for system identification. 
In the future we plan to compare the effectiveness of this 
methodology for system identification to that of other 
methods such as periodic driving [l3j and coupling a test 
system to a virtual model with tunable parameters [23| . 
The method we present need not be restricted to exam- 
ples such as these. In fact, the results are general and 
may be used to implement optimized control of any or 
all degrees of freedom in systems of ordinary differential 
equations. 
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t = t(p) with i(0) 

x i (t) = x i [t(p)]=x i (p), 



and t(l) 



x\t) 



(A.2) 
(A.3) 

(A.4) 



where a subscripted p indicates the partial deriva- 
tive with respect to the parameter p. By the scal- 
ing law for Dirac delta functions, f Q K5d [t{p) — r\dt = 

Jq KSD[t(p)— r]tpdp = J Q K5D(p—l)dp. Thus the func- 
tional then assumes the form 



SS = S I L{x\^-,t)t p + XK{x\^-,t)5 D {p- l)dp. 

tp tp 







(A.5) 



Then we execute the variation for each variable: 




OL 



dK 



Ox 



dK 



7717 *p + A 7717 <Mp ~ !) 



+ 



dL 



dK 
~dt 



-t p + X—6 D ( P -l) 



dx l 
Sx'l 

St 



dt p = 0. 



8x l 



(A.6) 



Next we evaluate the Dirac delta functions and use inte- 
gration by parts to eliminate 5x l and Sft. From Eq. (|A.6[) 
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we obtain: 



^=g{+(4f fai )L + f(s^* 



)=1 



J 



1 d f dL f 
dp y 

UdL 



p=l 



]-jfs( 1+ ^)«*- a 

(A.7) 



This gives rise to the following Euler-Lagrange equations 
for < p < 1: 

<9L d ( dL\ n ,. . 

^-*(^J=°' (A - 8) 
dL d / r <9L\ „ , A „, 

^flr-*( L+ ^J= ' (A - 9) 

for alH = 1, . . . , N. At the upper boundary, that is, for 
p = 1, we have 



dK 
~dT 



0, 



p 



dK T dL n 

x ^t +L + ^df p = °' 



air 



0, 



p 



t dL_ 

dx 1 p dx p 



(A.10) 
(A.ll) 
(A.12) 
(A.13) 



for all i = 1, . . . , JV. At the lower boundary, that is, 
for p = 0, we obtain for the variables where the initial 



conditions x l (0) and i l (0) are not fixed: 

<9L 



(A.14) 



Now we transform back to a parameter-free representa- 
tion with the following substitutions: 



x i (p)=x i [t(p)]=x i (t), 
4 = x l [t(p)]t p = x l (t)t p 



(A.15) 
(A.16) 



for all i = 1, . . . , N. Eqs. (|A~8]) and fO]) yield the same 
equation of motion: 







(A.17) 



dL d (dh 

dx l dt \ dx 1 
for all i = 1, . . . , N. At the upper boundary, that is, for 
t = t, we obtain from Eqs. (|A.10j) . (|A.11|) . (|A.12[) . and 
(IXTOl) . 



dK 



N r 



idL 
dx 1 



A 



^ = 0, 

dK dL 

dx 1 dx 1 



= 0, 



(A.18) 

(A.19) 
(A.20) 



for all i = 1, . . . , JV. At the lower boundary, that is, for 
t = t, for the variables where the initial conditions £*(()) 
and a; l (0) are not fixed, we obtain from Eq. (|A.14|) : 



dL 

dx* 



= 0, 



(A.21) 



Now we can use the Lagrange function given in Eq. ^ 
and substitute the variables x % in Eqs. (|A.17[) - (|A.21[) by 
x 1 ,. . .x d = x 1 ,. . . x d and x d+1 , . . . x 2d = F 1 , . . . F d , with 
N = 2d. 
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